#include "GBProcessing.h"

GBProcessing::GBProcessing()
{
    DBSeqNum=0;
}

GBProcessing::~GBProcessing()
{
    map<string, map<string,string> >::iterator it=Taxonomy.begin();
    while(it!=Taxonomy.end())
    {
        (*it).second.clear();
        it++;
    }
    Taxonomy.clear();
    TaxonAbundance.clear();
    SeqAbundance.clear();
}

string GBProcessing::get_Taxonomy(string SeqID, int tLevel)
{
    for(int i=tLevel;i>=0;i--)
        if(Taxonomy[SeqID].find(Level[i])!=Taxonomy[SeqID].end())
        {
             tLevel=i;
             break;
        }
    //if(Taxonomy[SeqID][Level[tLevel]].empty())
        //cout<<SeqID<<" has no taxonomy\n";
    return Taxonomy[SeqID][Level[tLevel]];
}

char* GBProcessing::string_to_char(string s)
{
    int len=strlen(s.c_str());
    char* temp=new char[len+1];
    memset(temp,0x00,sizeof(char)*(len+1));
    memcpy(temp,s.c_str(),len);
    return temp;
}

string GBProcessing::post_proc(string s)
{
    int len=s.size(), PrefixSpace=0;
    for(int i=0; i<s.size(); i++)
    {
        if(s[i]==' ')
        {
            len--;
            PrefixSpace++;
        }
        else
            break;
    }
    len+=(s[s.size()-1]=='.')? -1:0;
    s=s.substr(PrefixSpace, len);
    return s;
}

//********************* New Taxonomy Loading ********************
int GBProcessing::extract_taxonomy(string GBFileName)
{
    string GBLine,SeqID;
    map<string,string> SeqTaxon;
    vector<string> TaxonGBLines;
    fstream GBFile;
    cout<<"constructing taxonomy for "<<GBFileName<<endl;
    GBFile.open(GBFileName.c_str());
    if(GBFile.fail()){cout<<"cannot open GenBank File\n"; return 0;};
    while(!GBFile.eof())
    {
        getline(GBFile, GBLine);
        TaxonGBLines.clear();
        Tokenize(GBLine, TaxonGBLines, "\t");
        SeqID=TaxonGBLines[0];
        if(SeqID=="S002439377")
            for(int i=0;i<TaxonGBLines.size();i++)
                cout<<TaxonGBLines[i]<<"\n";
        int tTemp=0;
        if(TaxonGBLines.size()%3!=1){cout<<TaxonGBLines.size()<<"\t"<<SeqID<<" do not have 3*x members\n"; tTemp=1;}
        for(int i=1+tTemp;i<TaxonGBLines.size();i++)
        {
            //bool flag=false;
            //for(int j=0;j<9;j++)
               //if(TaxonGBLines[i+1]==Level[j]){flag=true; break;}
            //if(flag==true)
            //{
                SeqTaxon[TaxonGBLines[i+1]]=TaxonGBLines[i];
                i=i+2;
            //}
            
            /*if(flag==false)
            {
                int Current=i;
                do
                {
                    i++;
                    for(int j=0;j<9;j++)
                        if(TaxonGBLines[i+1]==Level[j]){flag=true; break;}                
                }while(flag==false);
                for(int j=Current;j<i;j++)
                {
                    SeqTaxon[TaxonGBLines[i+1]]+=TaxonGBLines[j];
                    if(j!=i-1)
                       SeqTaxon[TaxonGBLines[i+1]]+=" ";
                    i=i+2;
                }
            }*/
        }
        if(SeqID.size()!=0&&SeqTaxon.size()!=0)
        {
            DBSeqNum++;
            Taxonomy[SeqID]=SeqTaxon;
            SeqTaxon.clear();               
        }
    }
    GBFile.close();
    cout<<"Taxonomy Construction Complete\n";
    return 1;
}



/*
int GBProcessing::extract_taxonomy(string GBFileName)
{
    string GBLine,SeqID;
    map<string,string> SeqTaxon;
    fstream GBFile;
    cout<<"constructing taxonomy for "<<GBFileName<<endl;
    GBFile.open(GBFileName.c_str());
    if(GBFile.fail()){cout<<"cannot open GenBank File\n"; return 0;};
    while(!GBFile.eof())
    {
        getline(GBFile, GBLine);
        if(GBLine.substr(0,5)=="LOCUS")       //extract RDP ID
        {
            char* pch=NULL;
            char* temp=GBProcessing::string_to_char(GBLine);
            pch=strtok(temp," ");
            pch=strtok(NULL," "); //get RDP ID, which is the second component of LOCUS line
            SeqID.assign(pch);
            //cout<<SeqID<<"\n";
            delete[] temp;
        }
        if(GBLine.substr(0,10)=="  ORGANISM")  //extrac taxonomy
        {
            vector<string> TaxonGBLines;
            do
            {
                getline(GBFile, GBLine);
            }while(GBLine.substr(0,17)!="            Root;");
            TaxonGBLines.push_back(post_proc(GBLine));
            while(GBLine[GBLine.size()-1]!='.')  //PrevGBLine contains the genus-level taxonomy
            {
                getline(GBFile,GBLine);
                TaxonGBLines.push_back(post_proc(GBLine));
            }
            int nLines=TaxonGBLines.size();
            if(TaxonGBLines[nLines-1].size()<3)
            {
                TaxonGBLines[nLines-2]+=TaxonGBLines[nLines-1];
                TaxonGBLines.pop_back();           //if "." forms a single line, eliminate it
            }
            nLines=TaxonGBLines.size();
            //determine if unclassified:
            vector<string> Components;
            Tokenize(TaxonGBLines[nLines-1],Components,";.\r");
            size_t found1=Components[Components.size()-1].find("unclassified_");
            size_t found2=Components[Components.size()-1].find("Unclassified_");
            int CurrentLevel=0;
            if(found1!=string::npos||found2!=string::npos)  //unclassified
            {
                 if(nLines==3)
                 {
                     SeqTaxon[Level[8]]=post_proc(Components[0]); 
                     if(Components.size()!=1){cout<<"unclassified sequence "<<SeqID<<" has "<<nLines<<" lines but not only genus in last line\n";}
                     Components.clear();
                     Tokenize(TaxonGBLines[0],Components,";.\r");
                     for(int i=0;i<Components.size();i++) {SeqTaxon[Level[i]]=post_proc(Components[i]);}
                     Components.clear();
                     Tokenize(TaxonGBLines[1],Components,";.\r");
                     CurrentLevel=5;
                     for(int i=0;i<Components.size();i++)
                     {
                         if(Components.size()==2&&i==1)
                             CurrentLevel++; //no suborder
                         SeqTaxon[Level[i+CurrentLevel]]=post_proc(Components[i]);
                     }
                     Components.clear(); 
                 }
                 if(nLines==1)
                 {
                     Components.clear();
                     Tokenize(TaxonGBLines[0],Components,";.\r");
                     for(int i=0;i<Components.size();i++) {SeqTaxon[Level[i]]=post_proc(Components[i]);}
                     Components.clear();
                 }
                 if(nLines==2)
                 {
                     if(Components.size()==1) //only unclassified in 2nd line
                     {
                         Components.clear();
                         Tokenize(TaxonGBLines[0],Components,";.\r");
                         int tTemp=0;
                         for(int i=0;i<Components.size();i++) {if(i==4){tTemp=1;} SeqTaxon[Level[i+tTemp]]=post_proc(Components[i]); CurrentLevel=i+tTemp+1;}
                         if(CurrentLevel==4||CurrentLevel==6) CurrentLevel++;
                         Components.clear();
                         Tokenize(TaxonGBLines[1],Components,";.\r");
                         SeqTaxon[Level[CurrentLevel]]=post_proc(Components[0]);
                         Components.clear();
                     }
                     else
                     {
                         Components.clear();
                         Tokenize(TaxonGBLines[0],Components,";.\r");
                         int tTemp=0;
                         for(int i=0;i<Components.size();i++) {SeqTaxon[Level[i]]=post_proc(Components[i]); CurrentLevel=i+1;}
                         //if(SeqID=="S000021989"){cout<<Components.size()<<"\t"<<CurrentLevel<<"\n"; for(int i=0;i<Components.size();i++){cout<<Components[i]<<"1\n";}}
                         if(CurrentLevel==4) CurrentLevel++;
                         Components.clear();
                         Tokenize(TaxonGBLines[1],Components,";.\r");
                         for(int i=0;i<Components.size();i++) {if(Components.size()==2&&i==1){CurrentLevel++;} SeqTaxon[Level[i+CurrentLevel]]=post_proc(Components[i]);}
                         Components.clear();                     
                     }
                 }
            
            }
            else //classified
            {
                 
                 Components.clear();
                 for(int i=0;i<nLines;i++)
                     Tokenize(TaxonGBLines[i],Components,";.\r");
                 if(Components.size()>9)
                 {
                     cout<<SeqID<<"\t"<<Components.size()<<"\n";
                     for(int i=0;i<Components.size();i++)
                         cout<<Components[i]<<"\n";
                 }    
                 if(Components.size()<=7)
                 {
                     int tTemp=0;
                     for(int i=0;i<Components.size();i++)
                     {
                         if(i==Components.size()-1){ SeqTaxon[Level[8]]=post_proc(Components[i]); break;}
                         if(i==4||i==5)
                             tTemp++;
                         SeqTaxon[Level[i+tTemp]]=post_proc(Components[i]);
                     }
                 }
                 if(Components.size()>7)
                 {
                     int tTemp=0;
                     for(int i=0;i<Components.size();i++)
                     {
                         if(i==Components.size()-1){ SeqTaxon[Level[8]]=post_proc(Components[i]); break;}
                         if(i==5&&Components.size()==8)
                             tTemp++;
                         SeqTaxon[Level[i+tTemp]]=post_proc(Components[i]);
                     }
                 }
                 /*
                 if(nLines==2)
                 {
                     Components.clear();                 
                     Tokenize(TaxonGBLines[0],Components,";.\r");
                     int tTemp=0;                    
                     for(int i=0;i<Components.size();i++)
                     {
                         if(Components.size()>5)
                         {
                             if(i==4)
                                 tTemp=1;
                             if(i==5)
                                 tTemp=2;
                         }
                         SeqTaxon[Level[i+tTemp]]=post_proc(Components[i]);
                     }
                     Components.clear();
                     Tokenize(TaxonGBLines[1],Components,";.\r");
                     if(Components.size()>2) {CurrentLevel=5;} //if has order
                     else if(Components.size()==2){CurrentLevel=7;} //if no order but family
                     else if(Components.size()==1){CurrentLevel=8;}
                     for(int i=0;i<Components.size();i++)
                     {
                         if(Components.size()==3&&i==1)
                             CurrentLevel++; //no suborder
                         SeqTaxon[Level[i+CurrentLevel]]=post_proc(Components[i]);
                     } 
                     Components.clear();
                 }
                 if(nLines==1)
                 {
                     Components.clear();
                     Tokenize(TaxonGBLines[0],Components,";.\r");
                     for(int i=0;i<Components.size();i++) 
                     {
                         if(i==Components.size()-1)
                             SeqTaxon[Level[8]]=post_proc(Components[i]);
                         else
                             SeqTaxon[Level[i]]=post_proc(Components[i]);
                     }                     
                 }
                 if(nLines==3)
                 {
                     SeqTaxon[Level[8]]=post_proc(Components[0]); 
                     Components.clear();
                     Tokenize(TaxonGBLines[0],Components,";.\r");
                     for(int i=0;i<Components.size();i++) {SeqTaxon[Level[i]]=post_proc(Components[i]);}
                     Components.clear();
                     Tokenize(TaxonGBLines[1],Components," ;.\r");
                     CurrentLevel=5;
                     for(int i=0;i<Components.size();i++)
                     {
                         if(Components.size()==2&&i==1)
                             CurrentLevel++; //no suborder
                         SeqTaxon[Level[i+CurrentLevel]]=post_proc(Components[i]);
                     }
                     Components.clear();                   
                 }
                 
            }
        }
        if(SeqID.size()!=0&&SeqTaxon.size()!=0)
        {
            DBSeqNum++;
            Taxonomy[SeqID]=SeqTaxon;
            /*
            map<string, string>::iterator it=SeqTaxon.begin();
            cout<<SeqID<<"\t";
            while(it!=SeqTaxon.end())
            {
                cout<<(*it).first<<":"<<(*it).second<<"\t";
                it++;
            }
            cout<<"\n";
            
            SeqID.clear();
            SeqTaxon.clear();
        }
    }
    GBFile.close();
    return 1;
}
*/

int GBProcessing::load_taxonomy(string GBFileName)
{
    string TaxonomyFileName=GBFileName+".taxon";
    fstream TaxonFile;
    TaxonFile.open(TaxonomyFileName.c_str());
    if(TaxonFile.fail()){cout<<"no existing taxonomy file detected\n"; return 0;}
    string SeqID, TaxonLevel, SeqTaxon, Temp;
    TaxonFile>>DBSeqNum; //load database size
    getline(TaxonFile, Temp);
    while(!TaxonFile.eof())
    {
        TaxonFile>>SeqID>>TaxonLevel>>SeqTaxon;
        getline(TaxonFile,Temp);
        Taxonomy[SeqID][TaxonLevel]=SeqTaxon;
        map<string,int>::iterator it=TaxonAbundance.find(SeqTaxon);
        if(it==TaxonAbundance.end())
            TaxonAbundance[SeqTaxon]=1;
        else
            TaxonAbundance[SeqTaxon]++;
    }
    /*
    map<string, map<string, string> >::iterator it=Taxonomy.begin();
    while(it!=Taxonomy.end())
    {
        SeqAbundance[(*it).first]=TaxonAbundance[(*it).second];
        it++;
    }
    TaxonFile.close();
    */
    cout<<"Load Taxonomy Complete"<<endl;
    return 1;
}

int GBProcessing::output_taxonomy(string GBFileName)
{
    string TaxonomyFileName=GBFileName+".taxon";
    fstream TaxonFile;
    TaxonFile.open(TaxonomyFileName.c_str(),ios_base::out);
    if(TaxonFile.fail()){cout<<"cannot create taxonomy file"; return 0;}
    TaxonFile<<DBSeqNum<<"\n";
    map<string, map<string, string> >::iterator it=Taxonomy.begin();
    while(it!=Taxonomy.end())
    {
        map<string, string>::iterator it1=(*it).second.begin();
        while(it1!=(*it).second.end())
        {
            TaxonFile<<(*it).first<<"\t"<<(*it1).first<<"\t"<<(*it1).second<<"\n";
            it1++;
        }
        it++;
    }
    TaxonFile.close();
    return 1;
} 

/*
int GBProcessing::get_SeqAbundance(string SeqID)
{
   return SeqAbundance[SeqID]; 
}
*/

void GBProcessing::Iterator_reset(int Flag)
{
    if(Flag==1||Flag==0)
        GBIterator=SeqAbundance.begin();
    if(Flag==2)
        GBIterator=TaxonAbundance.begin();
    return;
}

bool GBProcessing::Iterator_end(int Flag)
{
    if(Flag==1||Flag==0)
    {
        if(GBIterator==SeqAbundance.end())
            return true;
        else
            return false;
    }
    if(Flag==2)
    {
        if(GBIterator==TaxonAbundance.end())
            return true;
        else
            return false;
    }
}

int GBProcessing::Iterator_size(int Flag)
{
    if(Flag==1||Flag==0)
    {
        return SeqAbundance.size();
    }
    if(Flag==2)
    {
        return TaxonAbundance.size();
    }
}

void GBProcessing::Tokenize(const string& str, vector<string>& tokens, const string& delimiters = " ")
{
    // Skip delimiters at beginning.
    string::size_type lastPos = str.find_first_not_of(delimiters, 0);
    // Find first "non-delimiter".
    string::size_type pos     = str.find_first_of(delimiters, lastPos);
    
    while (string::npos != pos || string::npos != lastPos)
    {
        // Found a token, add it to the vector.
        if(str.substr(lastPos, pos - lastPos)!=" ")
            tokens.push_back(str.substr(lastPos, pos - lastPos));
        // Skip delimiters.  Note the "not_of"
        lastPos = str.find_first_not_of(delimiters, pos);
        // Find next "non-delimiter"
        pos = str.find_first_of(delimiters, lastPos);
    }
    return;
}
